clear;
close all;
clc;
pause(0.2)
rad=0.15;
thre = 0.3;
% thre = 0.02;
% freqs= [300 600 1500 3000]; %old
% freqs= [850 1000];
% freqs= [850 900 1000];
freqs= 100:100:3000; %just a random frequency
fs = 44100;
Proto_Num = 5;
tdur=40e-3;
% out_table=f_P3_2_table;
% Proto_Num = 4.2;
% out_table=f_P4_2_table;
% Proto_Num = 5.2;
% out_table=f_P5_2_table;

if Proto_Num==3
    load 'C:\Users\gaiy\OneDrive - Saint Louis University\Sound_loc\rand_locations\P3_28-Jul-2025'
elseif Proto_Num==4
    load 'C:\Users\gaiy\OneDrive - Saint Louis University\Sound_loc\rand_locations\P4_25-Jul-2025'
elseif Proto_Num==5
    load 'C:\Users\gaiy\OneDrive - Saint Louis University\Sound_loc\rand_locations\P5_30-Jul-2025'
end
out_table=sort(all_locs')';
% ser0

% true_angle = [-80 -65 -55 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 55 65 80];
Elevation_number = 9;

% Angles = 1:3;
% Angles = 1:30;

conds = 'AB';
% sides = {'Left'}; % fake for Protocol = 3

PString = num2str(Proto_Num);
% f2=figure('position',[400 70 500 500]);
for iangle = 1:30
    % for iangle = 1:3
    cc = 0;
    f1=figure('position',[10 40 2100 900]);
    %     for icond = 1:1
    [recorded_sound,fs]=audioread(['C:\Users\gaiy\OneDrive - Saint Louis University\Sound_loc\3Dio recordings\P' PString '_Trial_' num2str(iangle) '.mp3']);
        recorded_sound=resample(recorded_sound(fs*1:end-1*fs,:),441,320);
        
    left_sound = recorded_sound(:,1);
    right_sound = recorded_sound(:,2);
    for ifreql = 1:length(freqs)
        freq = freqs(ifreql)    
        % tdur = 40/freq;
        [SL,SR,freq_pos,flag]=mySTFT(left_sound,right_sound,fs,thre,tdur);
        % normalization:
        [SL1,SR1]=Normalize(SL,SR);
        [Y Ind] = min(abs(freq_pos - freq));
        X_P = [real(SR1(Ind,flag>0))' imag(SR1(Ind,flag>0))'];
%             X_P = f_cleanup_Bnoise(X_P,fs,freq,tdur);

        cc=cc+1;
        figure(f1)
        subplot(5,6,cc)
        % for i=1:length(X_P)
        %     plot(X_P(i,1),X_P(i,2),'k.');hold on;
        % end
            plot(X_P(:,1),X_P(:,2),'k.');hold on;
        deg_sign=char(0176);

        out_center = f_mycenter(X_P,rad);
        %             sdf
        C = [out_center.C1; out_center.C2];
        % eval(['load Freq_centers' num2str(freq) ])
%         eval(['load Freq_centers_own' num2str(freq) ])
            eval(['load Model_centers\Freq_centers_sinu' num2str(freq) ])

        % for i=1:length(F_centers)
        %     s=plot(F_centers(i,1),F_centers(i,2),'go');hold on;set(s,'linewidth',2)
        % end
        ss=text(0.6,-0.8,[num2str(freq) ' Hz']);set(ss,'fontsize',14,'color','r')

        plot(C(1,1),C(1,2),'bo','Markersize',20,'LineWidth',2);hold on;
        plot(C(2,1),C(2,2),'mo','Markersize',20,'LineWidth',2);hold on;
        grid on
        set(gca,'Xtick',-1:0.5:1)
        set(gca,'Ytick',-1:0.5:1)
        axis([-1 1 -1 1])
        hold on;
        % plot(-1:1,[0,0,0],'k', 'LineWidth', 0.5);
        % plot([0,0,0],-1:1,'k','LineWidth',0.5)
        set(gca,'fontsize',10)
        xlabel('Re','fontsize',12)
        ylabel('Im','fontsize',12)
        hold on
        % now plot the line:
        % plot(x_spiral,y_spiral,'color',[0.5 0.5 0.5]) ; % nturns crossings, including end point
        axis square
        out = f_spiral_classification(C(1,:),x_spiral,y_spiral);
        % fine_angles = linspace(true_angle(1), true_angle(end), length(x_spiral));
        fine_angles = linspace(-90, 90, length(x_spiral));
        classified_angle(1)= round(fine_angles(out.class))

        out = f_spiral_classification(C(2,:),x_spiral,y_spiral);
        % fine_angles = linspace(true_angle(1), true_angle(end), length(x_spiral));
        fine_angles = linspace(-90,90, length(x_spiral));
        classified_angle(2)= round(fine_angles(out.class))
        %             record_class{icond}(ifreql,iangle)=classified_angle;
        if classified_angle(1) > classified_angle(2)
            classified_angle=fliplr(classified_angle)
            C = flipud(C)
        end
        % classified_angle=sort(classified_angle);
        title(['Classified:' num2str(classified_angle(1)) ' & ' num2str(classified_angle(2))])

        if ifreql == 1
            t = text(-3.5, 0.25, ['True:' num2str(out_table(iangle,1)) ' & ' num2str(out_table(iangle,2)) ]);
            set(t,'fontsize',9, 'FontWeight','bold')
        end
        record_class{ifreql,iangle}=[out_table(iangle,1:2) classified_angle];

        X_sp_freq(ifreql,:)=x_spiral;
        y_sp_freq(ifreql,:)=y_spiral;
            C_freq1(ifreql,:)=C(1,:);
            C_freq2(ifreql,:)=C(2,:);
     
        pause(0.1)
    end %ifreql
    %     end
    outMul = f_spiral_classification_MultiF(C_freq1,X_sp_freq,y_sp_freq)
    % fine_angles = linspace(true_angle(1), true_angle(end), length(x_spiral));
    fine_angles = linspace(-90,90, length(x_spiral));
    classified_MultiF1= fine_angles(outMul.class)

    outMul = f_spiral_classification_MultiF(C_freq2,X_sp_freq,y_sp_freq)
    % fine_angles = linspace(true_angle(1), true_angle(end), length(x_spiral));
    classified_MultiF2= fine_angles(outMul.class)
    record_classM{iangle}=[out_table(iangle,1:2) classified_MultiF1 classified_MultiF2];
    pause(1)
    pause
end

figure;
pc = 'bmgr';
for ifreql = 1:length(freqs)
    subplot(2,2,ifreql)
    curve = [];
    for iangle = 1:length(Angles)
    % for iangle = 1:3
        curve=[curve
            record_class{ifreql,iangle}(1) record_class{ifreql,iangle}(3)
            record_class{ifreql,iangle}(2) record_class{ifreql,iangle}(4)
            ];

    end
    [Y Ind ] =sort(curve(:,1));
    curve = curve(Ind,:)
    s=plot(curve(:,1),curve(:,2),['o' pc(ifreql)]);hold on;
    set(s,'linewidth',2)

    plot([-100 100],[-100 100],'k');
    % legend('300 Hz',' 600 Hz',' 1500 Hz',' 3000 Hz','Perfect','Location','northwest')
    % legend boxoff
    xlabel('True Angle (degrees)','fontsize',12)
    ylabel('classified Angle (degrees)','fontsize',12)
    title(['Dual Sources, ' num2str(freqs(ifreql)) ' Hz'])
    axis square
end


subplot(2,2,4)
curve = [];
for iangle = 1:length(Angles)
    % for iangle = 1:3
    curve=[curve
        record_classM{iangle}(1) record_classM{iangle}(3)
        record_classM{iangle}(2) record_classM{iangle}(4)
        ];
end
[Y Ind ] =sort(curve(:,1));
curve = curve(Ind,:)
s=plot(curve(:,1),curve(:,2),['or' ]);hold on;
set(s,'linewidth',2)
plot([-100 100],[-100 100],'k');
% legend('300 Hz',' 600 Hz',' 1500 Hz',' 3000 Hz','Perfect','Location','northwest')
% legend boxoff
xlabel('True Angle (degrees)','fontsize',12)
ylabel('classified Angle (degrees)','fontsize',12)
title(['Dual Sources, combining frequencies'])
axis square